Packages you will need:
library(RColorBrewer)
library(leaflet)
library(ggplot2)
library(mapview)
library(raster)
library(htmltools)
library(sf)
There are many ways to read spatial data into R
Here we get admin level 2 boundaries for Nepal from GADM (using raster package) and then convert the objcect to an sf class.
If you don’t have internet, you can read in shapefiles from your drive (as you might do in ArcGIS). R-packages also come with inbuilt spatial datasets- see this Mapview tutorial which uses an inbuilt dataset on breweries in a region of Germany
Nepal2 <- getData('GADM', country='NPL', level=2)
Nepal2 <- st_as_sf(Nepal2) # convert to an sf object (for later processing)
mapview(Nepal2)
# colour districts by name
mapview(Nepal2, zcol = "NAME_2")
Let’s generate some dummy data:
## generate random number for cases
Nepal2$cases <- sample(0:100, nrow(Nepal2), replace=T)
mapview(Nepal2, zcol = "cases")
It’s easy to define breaks using seq() function
mapview(Nepal2, zcol = "cases", at = seq(from = 1, to = 100, by =10), legend = TRUE)
Create a new field to show (made-up) rate
Nepal2$population <- sample(10000:30000, nrow(Nepal2), replace=T)
Nepal2$rate <- round(Nepal2$cases/Nepal2$population*10000,1)
mapview(Nepal2, zcol = "rate")
# generate some made up health facilities
facilities <- st_sample(x = Nepal2, size = 50, type = "random", exact = TRUE)
## plot Nepal districts and health facilities
mapview(Nepal2, zcol = "cases") + facilities
## change colour of facilities
mapview(Nepal2, zcol = "cases") + mapview(facilities, col.regions = "red")
# define your own palette to colour the areas
pal <- colorRampPalette(c("white", "brown1"))
mapview(Nepal2, zcol = "cases", at = seq(0, 100, 20), col.regions = pal(5), legend = TRUE)
Leaflet has many more options to customise the presentation of your map
The learning curve is steeper than for mapview (but there’s lots of online examples)
It’s ideal if you want to produce maps to present on a website
# Use leaflet to display a blank map with OpenStreetMap imagery
map <-leaflet() %>% ## using the pipe operator to chain operations together: here the blank map is chained to
addProviderTiles(providers$OpenStreetMap) ## the first leaflet option we are applying
map
Add Nepal districts to the OpenStreet base map
Map_Nepal<-map %>%
addPolygons(data=Nepal2)
Map_Nepal
colorBin() is a leaflet function that maps continuous numerica data onto a palette First make a gradient with as many colours as you like e.g. red-yellow-green
Then produce a palette mapping the values in your dataframe to this gradient
Or you can map the data to one of R’s pre-defined palettes (e.g. “Blues”)
My_gradient <- c("white", "brown1" ) ## brown1 is actually red
palTotal <- colorBin(My_gradient, domain = Nepal2$cases, bins = quantile(Nepal2$cases))
palRate <- colorBin("Blues", domain = Nepal2$rate, bins = quantile(Nepal2$rate))
# we will call these palettes in options when we make the map
This generates a simple, non-interactive map
Here we use different options to define the outline and fill style, and assign the layer to a group
There are numerous other options for you to explore here For example, try making the lines dashed and white instead of solid and black
Map_Nepal <- map %>% ## use the pipe operator to add elements
## to the OpenStreet base map
addPolygons(data = Nepal2, ## Use data from the spatial polygons dataframe Nepal2
color = 'black', ## line color
weight =2, ## line width
opacity = 1, ## line opacity
fillColor = ~palTotal(cases), ## fill polygons using the palette we generated
fillOpacity = 1, ## opaque fill
group = "Number of cases") ## allows this element to be linked e.g. to legend
Map_Nepal
When we start adding more layers to the map, we will need a control panel so the user can choose which layers to display
We will also need a layers panel that should be linked to the data and the control panel
Map_Nepal <- Map_Nepal %>% ## again we use the pipe operator
## to add complexity to our existing map object
## E.g. add the "Number of cases" group to the layers panel.
addLayersControl(overlayGroups = c("Number of cases")) %>%
## Currently this group only contains the polygons showing number of cases.
## Add a Legend
addLegend(pal = palTotal, ## use same palette as for polygons
values = Nepal2$cases, ## use same values
opacity = 1, ## self-explanatory
title = "Number of cases", ## "
position = "bottomleft", ## "
group = "Number of cases", ## links this legend with the polygons
labFormat = labelFormat(digits = 0)) ## define decimal places
Map_Nepal
We can add layers showing other fields just by changing the fillColor argument and assigning it to a new group
To make your maps look consistent, its usually best to keep the other options (e.g. line style) the same
Map_Nepal <- Map_Nepal %>%
addPolygons(data = Nepal2,
weight = 2,
opacity = 1,
color = 'black',
fillColor = ~palRate(rate), ## colour by rate
fillOpacity = 1,
group = "Rate") %>% ## new group
# Layers control
addLayersControl(overlayGroups = c("Number of cases", "Rate"), ## as before but with 2 groups
options = layersControlOptions(collapsed = FALSE))%>% ## this keeps the panel expanded
# Legend
addLegend(pal = palRate,
values = Nepal2$rate,
opacity = 1,
title = "Rate (per 10,000)",
position = "bottomleft",
group = "Rate",
labFormat = labelFormat(digits = 1))
Map_Nepal
This is a long chunk of code but it is just putting together what we have done above, and adding highlight and label options.
You could also look at this Example from leaflet
## define labels using the sprintf() function to format text
number_labels <- sprintf(
"<strong>%s</strong><br/>%d cases", # this is C-style formatting which means
Nepal2$NAME_2, Nepal2$cases # "print NAME_2 in bold font, return,
# print the number of cases in normal font"
) %>% lapply(htmltools::HTML)
rate_labels <- sprintf(
"<strong>%s</strong><br/>%g cases per 10,000",
Nepal2$NAME_2, round(Nepal2$rate, digits = 1)
) %>% lapply(htmltools::HTML)
## Add these labels to the map, along with highlighting options
## Number layer - inital options as above
Map_Nepal <- map %>%
addPolygons(data=Nepal2,
fillColor= ~palTotal(cases),
fillOpacity = 1,
color = 'black',
weight =2,
opacity = 1,
group = "Number of cases",
highlight = highlightOptions( # When user hovers, give the polygon:
weight = 4, # A thicker outline
color = "#666"), # A grey outline
label = number_labels, # The label we defined as number_labels
labelOptions = labelOptions(
textsize = "15px")) %>%
## rate
addPolygons(data=Nepal2,
fillColor= ~palRate(rate),
fillOpacity= 1,
color= 'black',
weight=2,
opacity= 1,
group = "Rate",
highlight = highlightOptions(
weight = 4,
color = "#666"),
label = rate_labels,
labelOptions = labelOptions(
textsize = "15px")) %>%
# Layers control
addLayersControl(overlayGroups = c("Number of cases", "Rate"),
options = layersControlOptions(collapsed = FALSE))%>%
# Legend
addLegend(pal = palTotal,
values = Nepal2$cases,
opacity = 1,
title = "Number of cases",
position = "bottomleft",
group = "Number of cases",
labFormat = labelFormat(digits = 0))%>%
addLegend(pal = palRate,
values = Nepal2$rate,
opacity = 0.7,
title = "Rate (per 10,000)",
position = "bottomleft",
group = "Rate",
labFormat = labelFormat(digits = 1))%>%
hideGroup("Rate") ## This means only the first group will be displayed initially
Map_Nepal
If you are making the map for a funder, a programme or any other external partner, you may need to add their logo
The url option is really useful if your maps are going to appear on a separate webpage
Map_Nepal <- Map_Nepal %>%
addLogo(img = "http://blogs.lshtm.ac.uk/inthenews/files/2017/12/LSHTM-logo-bw-1.jpg",
url = "https://www.lshtm.ac.uk/",
position = "topleft",
width = 120,
height = 60) %>% addMeasure()
Map_Nepal